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ABSTRACT 


Following  introduction,  tbe  paper  addresses  vertical  deflection 
estimation  methods  developed  by  Jordan  and  by  White  and  Goldstein, 
presents  the  development  of  two  independent  linear  least  squares 
estimation  methods  with  one  having  optimal  characteristics  as  to 
accuracy,  economy,  and  versatility,  and  arrives  at  significant 


conclusions. 


1.  INTRODUCTION 


Bell  Aerospace  completed  development  of  the  ball-bearing 
rotating  accelerometer  gravity  gradiometer  instrument  (6GI)  in 
1981.  During  that  year.  Bell  also  concluded  investigations 
concerning  gimbal  and  survey  vehicle  self  gradients.  GGI  tests  on 
board  a  ship  have  been  encouraging  and  have  yielded  data  for  the 
establishment  of  operational  self-noise  power  spectral  densities. 
Essentialy  during  the  period  1983-84,  Bell  accomplished  some  GGI 
modifications  for  airborne  and  terrestrial  applications.  A  variety 
of  corresponding  tests  will  be  conducted  during  1985  and  1986. 
Self-noise  power  spectral  densities  for  the  GGI  mounted  at  the 
umbrella  angle  (spin  axis  35®  from  horizontal)  have  been  approxi¬ 
mated  in  the  open  literature  by  4$  (f)  ■  “  Af +  6  vhere.  f  denotes 
frequency  in  cycles  per  second,  A  is  a  constant  relating  to  low 
frequency  or  red  noise,  and  B  is  a  constant  indici^ing  high 
frequency  or  white  noise.  White  (1980)  listed  A  *.2.10"®  E^'Hz  and 
16.10"®  E^*  Hz  and  B  =  81  E^*Hz"*  and  86  E^  Hz"^  for  configura¬ 


tions,,  wi„th  vertical  and  horUontal  spin  axis,  respectively  where  E 
is  Eotvos  unit  *  10"^  sec"^,  and  H  is  Hertz  =  1  sec"^.  As  to 


environmental  sensitivities,  linear  acceleration  influences  are  the 
most  severe,  particularly  at  multiples  of  the  nominal  rotation 
frequency  {]  =  0.25  Hz.  The  low-frequency  error  may  be  essentially 
interpreted  as  random  drift.  Non-spectral  data  processing  methods 
require  a  different  analytical  representation  of  along  track  and 
cross  track  self-noise  power  spectra  in  order  to  obtain  covariance 
functions  by  the  Fourier  cosine  transform.  Essentially  different 
moving  GGI  information  becomes  available  every  two  seconds. 
Gradient  filters  are  nominally  fourth  order  Butterworth  filters  to 
attenuate  in,  3n  and  other  noise  manifestations.  Signal  informa¬ 
tion  is  obtained  as  a  ten  second  moving  window  average.  Elimina¬ 
tion  of  systematic  drifts  and  biases  has  to  be  accomplished  at  the 
start  of  a  survey  mission.  In  this  respect,  initial  gravity  tensor 
information  would  be  valuable.  Other  preprocessing  includes  trans¬ 
formation  of  data  from  umbrella  to  local  level  frame,  residual 
compensations,  prefiltering,  removal  of  spikes  where  indicated,  and 
subtraction  of  normal  gravity  gradients.  Provided  that  sufficient 
regional  gravity  anomaly  and  vertical  deflection  data  are  avail¬ 
able,  regional  means  may  be  removed  and  variances  and  a  distance 
correlation  parameter  may  be  adjusted  at  the  outset.  Regionally 
improved  covariance  functions  should,  however,  only  apply  to  a 
gradiometer  data  generation  domain  sufficently  embedded  in  the 
region.  This  restriction  applies  also  if  gravity  anomalies  and 
vertical  deflections  are  estimated  by  means  of  standard  covariance 
functions  and  first  estimates  of  components  of  the  disturbance 


■ 


ft 


y  Codes 


d  !d/or 

wiJwCiai 


gravity  vector  are  utilized  for  covariance  parameter  improve¬ 
ments.  The  restriction  of  vertical  deflection  or  gravity  anomaly 
estimation  to  a  domain  somewhat  smaller  than  the  data  coverage  area 
is  further  based  on  estimation  analysis. 

Following  the  above  preliminaries,  this  paper  addresses 
successively  in  the  following  three  sections  Jordan's  spectral 
estimation  method  and  significant  conclusions,  the  two  stage 
estimation  method  developed  by  White  and  Goldstein  which  includes 
variable  template  averages  during  the  first  stage,  limitations 
associated  with  the  use  of  a^T/3z^-measurements,  and  two 
Independent  linear  least  squares  estimation  methods.  The  second 
has  optimal  characteristics  as  to  accuracy,  economy,  and 
versatility  both  to  vertical  deflection  and  gravity  anomaly 
determination.  As  a  result  of  this  research,  several  important 
conclusions  can  be  drawn. 

2.  JORDAN'S  SPECTRAL  ESTIMATION  METHOD  AND  GENERAL  FINDINGS. 


Jordan's  (1982)  spectral  estimation  method  is  confined  to 
gravity  anomaly  and  vertical  deflection  error  estimation.  Because 
of  the  assumption  of  steady  state  and  the  associated  availability 
of  measured  gradients  in  an  infinite  plane,  first  derivatives  of 
the  disturbing  potiential  T  and  not  differences  thereof  may  be 
estimated  by  measured  second  derivatives  of  T  all  of  which  contain 
a  noise  component.  Under  the  assumption  that  all  gravity-related 
quantities  are  stationary  random  variables  it  is,  for  example,  with 
the  functional  g  *  3T/3x  and  the  along  track  measurement  x  =  3^T/3x 
+  n  *  s  +  n, 


.1 

s  s 

gs  XX 


In  eq.  (1),  H  is  the  systems  function,  Sgs  is  the  siqnal  cross¬ 
spectrum,  and  Sxx  =  Sqj  +  S«n  -is  the  message  or  signal  plus  noise 
spectrum.  The  power  spectral  density  of  the  estimation  error  is 


Sgg(j<«>)  *  Sgg  -  H  Sgx  *  Sgg  -  HSg^  (2) 

The  variance  of  the  estimation  error  is  computed  by  integrating 
under  Sgg,  using  u  as  an  independent  variable.  In  this  respect, 
Jordan  excluded  wavelength  exceeding  500  Km.  In  the  case  of  two 
independent  gravity  gradients,  (jw)  and  H2  (Ju)  would  have  to  be 
determined  from 

V.  <“>  <3> 


Snv.  (J“)  ■  S.  (“)  ”1  (J“)  ♦  5 _ N,  (j») 


In  order  to  solve  eqs.  (3)  and  (4)  for  Hj  and  Ho,  noise  cross¬ 
spectra  would  be  required  as  well;  however,  this  tends  to  be 
neglected  for  simplicity.  Jordan  established  the  signal  and  cross¬ 
signal  spectra  from  a  three-dimensional  algebraic  gravity  model 
(STAG).  The  STAG  model  was  fitted  to  a  third-order  Markov  model 
and  covariance  functions  involving  T  and  first  and  second 
derivatives  of  T  were  established  by  Jordan,  Moonan  and  Weiss 
(1981).  Because  the  STAG  model  permits  the  use  in  three 
dimensions,  Jordan  estimated  errors  at  ground  level.  In  accordance 
with  the  above  analysis,  a  matrix  of  dimension  5  N  X  5  N  has  to  be 
inverted  if  five  independent  gradient  measurements  and  N  parallel 
tracks  are  employed  for  error  estimation.  Jordan  used  A  *  17.10"® 
Hz  and  B  values  from  30  to  630  Hz"^  for  the  computation  of 
standard  errors  under  inclusion  of  analytical  downward 
continuation.  Table  1  shows  single  track  gravity  component 
standard  errors  for  a  track  length  of  500  Km,  an  aircraft  speed  of 
300  Km  hr"^,  an  aircraft  altitude  of  600  m,  and  a  self-noise 
characterized  by  B  =  160  E‘  Hz"*. 


along-track 

cross-track 

vertical 


Standard  Errors  (mgal ) 
1.9 
2.2 
2.0 


TABLE  1.  SINGLE  TRACK  STANDARD  ERRORS 


The  relatively  high  single  track  standard  error  of  about  2  mgal 
or  0.4  arc  sec  is  due  to  the  wavelength  restriction  of  500  Km.  For 
this  reason,  gradiometer  measurements  generated  along  parallel  and 
cross-tracks,  i.e.,  a  greater  sample  size,  is  required  to  achieve 
smaller  standard  errors.  To  assure  strong  signal  correlations,  the 
track  spacing  should  be  small  and  economical.  Utilization  of 
cross-track  data  for  the  same  number  of  parallel  and  cross-tracks 
results  in  an  error  reduction  by  about  (/z)"  .  Table  2  shows 
vertical  deflection  and  gravity  anomaly  standard  errors  due  to  the 
utilization  of  five  north-south  tracks  and  five  east-west  tracks 
and  a  track  spacing  of  5  Km  under  otherwise  unchanged  conditions. 


north  vertical  deflection 
east  vertical  deflection 
gravity  anomaly 


Standard  Errors 
0.12  arcsec 
0.12  arcsec 
0.63  m  gal 


TABLE  2.  MULTIPLE  TRACK  STANDARD  ERRORS 
(5  PARALLEL  TRACKS,  5  PARALLEL  CROSS-TRACKS) 


Smaller  variations  in  self-noise  (B),  aircraft  speed,  aircraft 
altitude,  and  track  spacing  from  the  parameters  used  do  not 
significantly  affect  the  standard  errors  shown  in  Table  2.  It 
must,  however,  be  emphasized  that  these  standard  errors  are  not 
uniformly  valid  in  the  survey  area  and  increase  towards  the 
boundary.  The  necessity  of  using  at  least  two  vertical  deflections 


and  one  gravity  anomaly  In  a  limited  domain  has  also  an  Impact  on 
nonuniform  error  estimation. 

3.  ESTIMATION  METHOD  OF  WHITE  AND  GOLDSTEIN 


A  square  survey  area  containing  61  flight  paths  parallel  to  one 
s.ac  with  a  path  length  of  300  Km  and  a  spacing  of  5  Km  and  61 
orthogonal  flight  paths  with  equal  length  and  spacing  provides  for 
a  total  track  length  of  36,600  Km.  If  five  Independent 
preprocessed  gradient  measurements  are  available  at  1  Km  Intervals, 
the  measurement  vector  would  comprise  183,000  numbers,  requiring 
the  Inversion  of  a  prohibitively  large  positive  definite  matrix. 
For  this  reason.  White  and  Goldstein  (1984)  proposed  the 
utilization  of  aggregated  gradients  obtained  by  averaging 
measurements  along  peripheries  of  squares  centered  above  the  point 
of  estimation,  evident  from  Figure  1. 


CENTtR  POINT 


FIGURE  1.  DATA  TEMPLATE  WINDOWS  FOR  DATE  AVERAGING 

As  an  Illustration,  the  authors  define  a  truth  vector  under 
restriction  to  T„  -  gravity  gradlometer  measurements  and  12 
template  zones  as:  _  _ 


In  eq.  (5),  r  designates  position  vector,  6  is  the  gravity 

disturbance  vector"  consisting  of  three  component’s,  and  6  is  the 

local  mean  disturbance  vector. 

The  residual  gravity  disturbance  vector  may  then  be  estimated 
in  the  form 

d(r)  -  6(r)  -  «  ^(r)  =  B  x  (6) 

where  the  3  x  n  output  matrix  B  selects  those  linear  combinations 

of  the  elements  of  x  that  are  being  estimated. 

A  solution  of  eq.  (6)  can  be  achieved  by  recursive  Kalman 
estimation  algorithm  under  limitations  associated  therewith.  In 
this  respect,  the  authors  entertain  the  possibility  of  reducing  red 
noise  covariances  by  utilization  of  gradiometer  measurements  at 
flight  traverse  crossing  points  without  an  analysis  to  achieve  such 
reduction.  The  quadratic  template  method  requires  a  modification 
in  the  neighborhood  of  the  survey  boundary.  This  restriction 
applies  to  some  extent  to  any  estimation  method. 

Although  eq.  (5),  augmented  by  additional  averaged  gradiometer 
measurements,  represents  a  convenient  truth  vector,  the  aggregate 
gradient  estimators  are  neither  best,  complete,  or  suboptiomal 
estimators.  As  discussed  in  the  following  section,  highly 
accurate,  economic,  and  versatile  determinations  of  d  (r)  under 
consideration  of  measurements  of  3^T/ax^,  3^T/3xay,  3^T/§y‘r  3^T/»x3z, 
and  a2T/3xay  require  the  utilization  of  at  least  one  measured 
gravity  disturbance  vector. 


4.  DEVELOPMENT  OF  TWO  INDEPENDENT  TWO-STAGE  LINEAR  LEAST-SOUARES 
ESTIMATION  METHODS 


The  utilization  of  a^T/3z^  -  measurements  without  knowledge  of 
the  components  of  the  disturbance  vector  on  the  boundary  does  not 
contribute  much  to  gravity  vector  component  determination  with 
desired  rms  errors  of  about  0.2  arcsec.  As  shown  by  Baussus  von 
Luetzow  (1977),  vertical  deflections  C  and  n  can  be  found  as  the 
solution  of  a  Poisson  equation  in  planar  coordinates  from  measure¬ 
ments  in  the  interior  of  a  region,  preferably  a  rectangle  or  a 
circle,  and  boundary  deflections.  In  the  case  of  a  quadrangle  with 
a  side  of  300  Km,  the  omission  of  boundary  values  leads  to  a 
computation  error  of  approximately  1.3  arcsec  for  the  center 
vertical  deflection  component.  For  an  extended  side  of  1,000  Km, 
the  error  would  amount  to  about  0.7  arcsec.  The  errors  tend  to 
Increase  towards  the  boundary.  In  ^the  context  of  a  parallel  track 
estimation  scheme,  the  emphasis,  therefore,  has  to  be  on  the 
determination  of  3T/ax  from  3^T/3x^  -  and  3^T/3x3y  measurements  and 
on  a  similar  determination  of  3T/3y  and  3T/3z.  This  requires  the 
use  of  at  least  one  set  of  measured  first  order  derivatives  of  T. 
In  this  respect,  it  should  be  emphasized  that  the  use  of  a  limited 
number  of  astrogeodetic  vertical  deflections  with  rms  errors  of 
0.10  arcsec  and  of  gravity  anomalies  is  not  uneconomic  and  is 
expected  to  Improve  estimations  in  regions  close  to  the  survey 
boundary.  A  utilization  of  more  than  the  above  gradiometer 


»  _  •  -  •  te  '  ^  »  »  m 


measurements  for  the  estimation  of  only  one  first  order  derivative 
of  T  does  not  significantly  reduce  the  estimation  error  and  Is 
uneconomic.  It  may,  however,  be  possible  to  consider  3^T/az^- 
measurements  for  an  Independent  computation  of  C  and  n  after  their 
initial  estimation  on  the  boundary.  A  statistical  estimation 

of  f  and  7  would  only  slightly  reduce  the  above  errors  8*T/3x*  - 
and  3*T/ax3y  measurements  along  both  “horizontal"  and  '"vertical" 
traverses  are  to  be  used  for  estimation;  and  data  relating  to  about 
seven  tracks  would  have  to  be  employed  to  reduce  the  estimation 
error,  as  Indicated  In  Figure  2. 


Measurements  of 


Measurements  of 


i  9  /C  //  /2  /J  J¥ 

FIGURE  2.  TRACK  SYSTEM  FOR  ESTIMATION  OF  |I  AT  Pg 
FROM  -  and  -  MEASUREMENTS 

In  a  two>stage  optimization,  with  gradlometer  measurements 
suitably  arranged  over  a  grid  length  of  5  Km  to  avoid  formidable 
matrix  inversions,  the  14  measurement  matrices  may  be  written  In 
the  compact  form 


with  1,  *  1,  2,  7;  I2  »  8,  9,  ...,  14;  K  “  1.  2,  60  for 

a  survey  guadrangle  with  a  side  of  300  Km.  The  symbol  *  indicates 
averaged  gradient. 

Fourteen  first-stage  solutions  may  be  obtained  from  covariance 
-equations  of  the  structure 

,b, 

where  X  |  ^  refers  to  the  first  2  matrix  elements  of  eq.  (7),  is 
the  matrix  of  120  regression  coefficients  for  a  specific  i,  the  bar 
symbolizes  the  covariance  operator,  and  T  indicates  the  transpose 
operator. 

The  second-stage  solution  can  then  be  found  from  the  regression 
equation 


■  biAj  +  b2A2  +  ...  +  bi4Am  =  BAi  '  ' 

Because  of  a  track  length  of  300  km  instead  of  the  track  length  of 
500  km  used  by  Jordan,  the  standard  errors  listed  in  Table  2  will 
be  exceeded.  For  some  points,  the  estimation  might  be  performed 
under  utilization  of  quasi -diagonal  parallel  lines  longer  than 
measurement  tracks.  This  involves,  however,  a  re-arrangement  of 
suitably  averaged  gradiometer  measurements.  The  use  of  an  ergodic 
estimation  procedure  including  the  use  of  covariance  functions  for 
long  distances  further  contributes  to  an  rms  error  increase  to 
0.20-0.25  arcsec.  The  corresponding  error  in  the  boundary  region 
with  width  15  Km  will  be  somewhat  greater. 

In  view  of  the  above  and  due  to  the  fact  that  initial  gravity 
anomaly  errors  are  small,  the  following  alternative  estimation 
method  is  presented  which  would  simultaneously  offer  estimates 
essentially  independent  of  the  Wiener-Kolmogorov  optimization 
characterized  by  eqs.  (8)  and  (9). 

The  alternative  two-stage  estimation  method  presupposes  the 
availability  of  highly  accurate  vertical  deflection  components  with 
rms  errors  of  0.1  arcsec  and  of  a  97/3^  -component  with  an  rms 
error  of  0.1  mgal  at  the  center  of  the  survey  square  with  tracks 
300  Km  long.  A  straight  line  is  drawn  between  the  center  or 
initial  point  Pq  and  the  estimation  point  Pg,  for  convenience  to 
coincide  with  a  track  intersection  point.  As  evident  from  Figure 
3,  the  connecting  line  intersects  with  at  most  30  grid  squares. 


FIGURE  3.  INTERSECTED  SQUARES  CONTAINING 
8  AVERAGE  MEASUREMENTS  PER  SQUARE. 

The  measurement  vector  for  a  square  ,  1  dent  1*f led  by  the 
superscript  1  rtay  be  written  as 


In  eq.  (10),  the  symbol  *  Indicates  an  average  measurement  over 
a  square  grid  length,  and  o,u,1,  and  r  upper,  lower,  left,  and 
right  sides  of  a  square  grid,  respectively. 

There  are  1  measurement  vectors  relating  to  the  line  P-  P.,  and 
there  are  71  measurement  vectors  If  six  parallel  lines  »nth  a 
separation  of  5  Km  are  drawn,  also  evident  from  Figure  3.  Except 
for  the  boundary  region,  this  Is  possible. 


As  to  first-stage  estimation,  the  regression  equation 

i  *  n,  k  *  * 

(T-T)=A=J  am  (11) 

X  X  »  »  ik  ik 

e  0  -  i  «  1,  k  *  1 

holds  with  maximally  240  coefficients  to  be  determined.  For 
greater  accuracy,  Pq  Pg  may  be  extended  to  Include  6  or  more 
additional  measurement  vectors. 

Equation  (11)  is  applied  with  respect  to  the  6  lines  parallel 
to  Pp  Pg,  so  that  altogether  the  solutions  (Txg  -  Txp)  u  *  An  are 
obtained,  where  u  =  1,  2,  •**  7  except  for  the  survey  boundary 
region. 

The  second-stage  estimation  yields  the  final  estimate 


(Txg-  Txq)  =  biAi  +  b2A2  +  _  +  b7&7  (12) 

where  the  regression  coefficients  bjj  are  computed  in  the  usual 
manner.  The  estimation  by  means  of  eqs.  (11)  and  (12)  is  optimal 
as  to  accuracy  and  efficiency,  and  the  error  associated  with  the 
second-stage  estimation  of  Tx  is  expected  not  to  exceed  the  rms 
error  of  0.18  arcsec  with  the  exception  of  the  boundary  region. 
Further,  this  estimation  is  relatively  less  dependent  on  long 
distance  covariances. 

To  meet  a  requirement  of  0.18  arcsec  rms  for  the  boundary 
region,  highly  accurate  disturbance  vector  components  would  be 
required  at  the  survey  area  ctyners,  or  the  survey  area  would  have 
to  be  enlarged  to  315  X  315  Km^. 

The  second  two-stage  estimation  method  developed  would  also 
facilitate  the  computation  of  survey  area  mean  values  3T/3x,  3T/3y, 
aT/a  and  the  improvement  of  a  priori  covariance  functions. 
Gravity  anomaly  coverage  in  9  adjacent  survey  areas  and  utilization 
of  a  spherical  harmonics  representation  of  Ag,  resulting  in  a  modi¬ 
fied  Stokes 'function,  would  allow  an  independent  determination  of  Kq 
and  in  the  central  survey  area. 

As  already  stressed  by  White  and  Goldstein,  a  considerable 
effort  has  to  be  made  as  to  data  Identification  and  organization. 
The  computational  effort  involved  requires  the  use  of  standardized 
survey  areas  and  track  separations.  Regression  coefficients  need 
only  be  computed  once  under  exploitation  of  symmetries. 

Immediate  analytical  downward  continuation  has  not  been 
considered  since  the  estimation  process  is  simpler  in  a  horizontal 
plane,  and  is  thus  not  dependent  on  the  Sperry  three-dimensional 
algebraic  gravity  model  (STAG).  Downward  continuation  can  be 
linearly  applied  by  means  of  representative  vertical  gradients 
determined  from  gradiometer  measurements  for  low  altitudes. 


""vn  ii<;ton$. 

Jordan's  analyses  are  basic  as  to  error  estimation  of  gravity 
vector  components  exclusively  from  measurements  of  gravity 
a.j^.ciiLs  under  consideration  of  important  parameters.  Error 
estimates  become*.— Iwwever^  greater  with  a  decreasing  quadratic 
survey  area, for  which  the  spectral  estimation  approach  is  not 
suitable.  In  \^ition,  the  ergodic  estimation  structure  employed 
is  less  promisin^'for^.thfi  estimation  of  point  vector  components 
than  for  the  estimation  of  vector  cofnponefttL^jJferences  because  of 
the  existence  of  an  underlying  initial  value  pr^bblem^.  j,White  and 
Goldstein  attempt  to  use  weakly  and  strongly  averaged  aggregates  of 
all  gradiometer  measurements,  covering  the  whole  survey  area,  for 
the  estimation  of^,  point  vector  components,  starting  with  the 
consideration  of  ■^measurements.  ,  As  to  the  utilization  of  other 
tensor  measurements  and  estimation  of  vector  components  at  points 
outside  an  interior  survey  area  region,  they  would  have  to  address 
the  problem  of  quasi-optimal  data  aggregation.  The  utilization  of 
J^z  -measurements  is  insufficient  for  an  accurate  estimation  of 
vertical  deflections.  A  two-stage  estimation  process  similar  to 
Jordan's  spectral  approach  but  applicable  to  a  track  length  of  300 
Km  and  requiring  a  relatively  simple  data  averaging  effort  has  been 
outlined.  -^lA  developed  two-stage  estimation  of  gravity  vector 
component  differences,  to  be  augmented  by  a  highly  accurate  gravity 
vector  in  the  center  of  the  survey  area,  is  most  promising  from  the 
standpoint  of  accuracy,  economy,  and  versatility.  Four  additional 
gravity  vector  components  at  the  corners  of  the  quadratic  survey 
area  would  alleviate  the  problem  of  accuracy  deterioration  close  to 
the  survey  area  boundary.  ^The  second  method  solves  the  problem  of 


optimal  data  aggregation'^and'  provides  estimates  essentially 
ioHononcient  of  those  obtained  by  the  first  method.  It  can  provide 
an  estimate  of  gravity  component  mean  values  relating  to  the  survey 
area  and  a  possible  improvement  of  a  priori  covariance  functions. 
Under  certain  conditions,  it  permits  the  computation  of  the 
vertical  deflection  components  -at  the  survey  area  center.  Both 
estimation  methods  do  not  use  spatial  covariance  functions,  a 
useful  simplification,  since  representative  analytical  downward 
continuation  of  gravity  vector  components  can  be  accomplished  by 
means  of  gradiometer  measurements  without  appreciable  error  for 
altitudes  of  the  order  600  m. 


aggregation 
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